Corona_Analysis
Steven Smith, PhD
3/18/2020
The 2019-2020 Coronavirus Pandemic Analysis
Contact: Smith Research
BACKGROUND & APPROACH
I wanted to track and trend the coronavirus outbreak on my own curiosity. There are some interesting questions that may fall out of this, as it is a very historic moment, including scientifically and analytically (we have a large amount of data being shared across the globe, analyzed in real-time). The world has come to a halt because of it.
This analysis attempts to answer the following questions (more to come):
- What does the trend of the pandemic look like to date?
- What are future case predictions based on historical model?
- What interesting quirks or patterns emerge?
ASSUMPTIONS & LIMITATIONS: * This data is limited by the source. I realized early on that depending on source there were conflicting # of cases. Originally I was using JHU data… but this was always ‘ahead’ of the Our World In Data. I noticed that JHU’s website was buggy- you clicked on the U.S. stats but it didn’t reflect the U.S.. So I changed data sources to be more consistent with what is presented in the media (and Our World In Data has more extensive plots I can compare my own to). An interesting aside might be why the discrepancy? Was I missing something?
* Defintiions are important as is the idea that multiple varibales accumulate in things like total cases (more testing for example).
SOURCE RAW DATA: * https://ourworldindata.org/coronavirus
* https://github.com/CSSEGISandData/COVID-19/
*
INPUT DATA LOCATION: github (https://github.com/sbs87/coronavirus/tree/master/data)
OUTPUT DATA LOCATIOn: github (https://github.com/sbs87/coronavirus/tree/master/results)
TIMESTAMP
Start: ##—— Sat May 16 11:05:31 2020 ——##
PRE-ANALYSIS
The following sections are outside the scope of the ‘analysis’ but are still needed to prepare everything
UPSTREAM PROCESSING/ANALYSIS
- Google Mobility Scraping, script available at get_google_mobility.py
# Mobility data has to be extracted from Google PDF reports using a web scraping script (python , written by Peter Simone, https://github.com/petersim1/MIT_COVID19)
# See get_google_mobility.py for local script
python3 get_google_mobility.py
# writes csv file of mobility data as "mobility.csv"
SET UP ENVIORNMENT
Load libraries and set global variables
# timestamp start
timestamp()
## ##------ Sat May 16 11:05:31 2020 ------##
# clear previous enviornment
rm(list = ls())
##------------------------------------------
## LIBRARIES
##------------------------------------------
library(plyr)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0 ✓ purrr 0.3.3
## ✓ tibble 3.0.0 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange() masks plyr::arrange()
## x purrr::compact() masks plyr::compact()
## x dplyr::count() masks plyr::count()
## x dplyr::failwith() masks plyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::id() masks plyr::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::mutate() masks plyr::mutate()
## x dplyr::rename() masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(ggplot2)
library(reshape2)
##
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
library(plot.utils)
library(utils)
library(knitr)
##------------------------------------------
##------------------------------------------
# GLOBAL VARIABLES
##------------------------------------------
user_name <- Sys.info()["user"]
working_dir <- paste0("/Users/", user_name, "/Projects/coronavirus/") # don't forget trailing /
results_dir <- paste0(working_dir, "results/") # assumes diretory exists
results_dir_custom <- paste0(results_dir, "custom/") # assumes diretory exists
Corona_Cases.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
Corona_Cases.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv"
Corona_Deaths.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv"
Corona_Deaths.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"
Corona_Cases.fn <- paste0(working_dir, "data/", basename(Corona_Cases.source_url))
Corona_Cases.US.fn <- paste0(working_dir, "data/", basename(Corona_Cases.US.source_url))
Corona_Deaths.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.source_url))
Corona_Deaths.US.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.US.source_url))
default_theme <- theme_bw() + theme(text = element_text(size = 14)) # fix this
##------------------------------------------FUNCTIONS
List of functions
| function_name | description |
|---|---|
| prediction_model | outputs case estumate for given log-linear moder parameters slope and intercept |
| make_long | converts input data to long format (specialized cases) |
| name_overlaps | outputs the column names intersection and set diffs of two data frame |
| find_linear_index | finds the first date at which linearaity occurs |
##------------------------------------------
## FUNCTION: prediction_model
##------------------------------------------
## --- //// ----
# Takes days vs log10 (case) linear model parameters and a set of days since 100 cases and outputs a dataframe with total number of predicted cases for those days
## --- //// ----
prediction_model<-function(m=1,b=0,days=1){
total_cases<-m*days+b
total_cases.log<-log(total_cases,10)
prediction<-data.frame(days=days,Total_confirmed_cases_perstate=total_cases)
return(prediction)
}
##------------------------------------------
##------------------------------------------
## FUNCTION: make_long
##------------------------------------------
## --- //// ----
# Takes wide-format case data and converts into long format, using date and total cases as variable/values. Also enforces standardization/assumes data struture naming by using fixed variable name, value name, id.vars,
## --- //// ----
make_long<-function(data_in,variable.name = "Date",
value.name = "Total_confirmed_cases",
id.vars=c("case_type","Province.State","Country.Region","Lat","Long","City","Population")){
long_data<-melt(data_in,
id.vars = id.vars,
variable.name=variable.name,
value.name=value.name)
return(long_data)
}
##------------------------------------------
## THIS WILL BE IN UTILS AT SOME POINT
name_overlaps<-function(df1,df2){
i<-intersect(names(df1),
names(df2))
sd1<-setdiff(names(df1),
names(df2))
sd2<-setdiff(names(df2),names(df1))
cat("intersection:\n",paste(i,"\n"))
cat("in df1 but not df2:\n",paste(sd1,"\n"))
cat("in df2 but not df1:\n",paste(sd2,"\n"))
return(list("int"=i,"sd_1_2"=sd1,"sd_2_1"=sd2))
}
##------------------------------------------
##------------------------------------------
## FUNCTION: find_linear_index
##------------------------------------------
## --- //// ----
# Find date at which total case data is linear (for a given data frame)
## --- //// ----
find_linear_index<-function(tmp,running_avg=5){
tmp$Total_confirmed_cases_perstate.log<-log(tmp$Total_confirmed_cases_perstate,2)
derivative<-data.frame(matrix(nrow = nrow(tmp),ncol = 4))
names(derivative)<-c("m.time","mm.time","cumsum","date")
# First derivative
for(t in 2:nrow(tmp)){
slope.t<- tmp[t,"Total_confirmed_cases_perstate.log"]- tmp[t-1,"Total_confirmed_cases_perstate.log"]
derivative[t,"m.time"]<-slope.t
derivative[t,"date"]<-as.Date(tmp[t,"Date"])
}
# Second derivative
for(t in 2:nrow(derivative)){
slope.t<- derivative[t,"m.time"]- derivative[t-1,"m.time"]
derivative[t,"mm.time"]<-slope.t
}
#Compute running sum of second derivative (window = 5). Choose point at which within 0.2
for(t in running_avg:nrow(derivative)){
slope.t<- sum(abs(derivative[t:(t-4),"mm.time"])<0.2,na.rm = T)
derivative[t,"cumsum"]<-slope.t
}
#Find date -5 from the stablility point
linear_begin<-min(derivative[!is.na(derivative$cumsum) & derivative$cumsum==running_avg,"date"])-running_avg
return(linear_begin)
}READ IN DATA
- total number of cases. current source: https://github.com/CSSEGISandData (precvious source https://ourworldindata.org/coronavirus)
# Q: do we want to archive previous versions? Maybe an auto git mv?
##------------------------------------------
## Download and read in latest data from github
##------------------------------------------
download.file(Corona_Cases.source_url, destfile = Corona_Cases.fn)
Corona_Totals.raw <- read.csv(Corona_Cases.fn, header = T, stringsAsFactors = F)
download.file(Corona_Cases.US.source_url, destfile = Corona_Cases.US.fn)
Corona_Totals.US.raw <- read.csv(Corona_Cases.US.fn, header = T, stringsAsFactors = F)
download.file(Corona_Deaths.source_url, destfile = Corona_Deaths.fn)
Corona_Deaths.raw <- read.csv(Corona_Deaths.fn, header = T, stringsAsFactors = F)
download.file(Corona_Deaths.US.source_url, destfile = Corona_Deaths.US.fn)
Corona_Deaths.US.raw <- read.csv(Corona_Deaths.US.fn, header = T, stringsAsFactors = F)
# latest date on all data:
paste("US deaths:", names(Corona_Deaths.US.raw)[ncol(Corona_Deaths.US.raw)])## [1] "US deaths: X5.15.20"paste("US total:", names(Corona_Totals.US.raw)[ncol(Corona_Totals.US.raw)])## [1] "US total: X5.15.20"paste("World deaths:", names(Corona_Deaths.raw)[ncol(Corona_Deaths.raw)])## [1] "World deaths: X5.15.20"paste("World total:", names(Corona_Totals.raw)[ncol(Corona_Totals.raw)])## [1] "World total: X5.15.20"PROCESS DATA
- Convert to long format
- Fix date formatting/convert to numeric date
- Log10 transform total # cases
##------------------------------------------
## Combine death and total data frames
##------------------------------------------
Corona_Totals.raw$case_type<-"total"
Corona_Totals.US.raw$case_type<-"total"
Corona_Deaths.raw$case_type<-"death"
Corona_Deaths.US.raw$case_type<-"death"
# for some reason, Population listed in US death file but not for other data... Weird. When combining, all datasets will have this column, but US deaths is the only useful one.
Corona_Totals.US.raw$Population<-"NA"
Corona_Totals.raw$Population<-"NA"
Corona_Deaths.raw$Population<-"NA"
Corona_Cases.raw<-rbind(Corona_Totals.raw,Corona_Deaths.raw)
Corona_Cases.US.raw<-rbind(Corona_Totals.US.raw,Corona_Deaths.US.raw)
#TODO: custom utils- setdiff, intersect names... option to output in merging too
##------------------------------------------
# prepare raw datasets for eventual combining
##------------------------------------------
Corona_Cases.raw$City<-"NA" # US-level data has Cities
Corona_Cases.US.raw$Country_Region<-"US_state" # To differentiate from World-level stats
Corona_Cases.US.raw<-plyr::rename(Corona_Cases.US.raw,c("Province_State"="Province.State",
"Country_Region"="Country.Region",
"Long_"="Long",
"Admin2"="City"))
##------------------------------------------
## Convert to long format
##------------------------------------------
#JHU has a gross file format. It's in wide format with each column is the date in MM/DD/YY. So read this in as raw data but trasnform it to be better suited for analysis
# Furthermore, the World and US level data is formatted differently, containing different columns, etc. Recitfy this and combine the world-level stats with U.S. level stats.
Corona_Cases.long<-rbind(make_long(select(Corona_Cases.US.raw,-c(UID,iso2,iso3,code3,FIPS,Combined_Key))),
make_long(Corona_Cases.raw))
##------------------------------------------
## Fix date formatting, convert to numeric date
##------------------------------------------
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "^X",replacement = "0") # leading 0 read in as X
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "20$",replacement = "2020") # ends in .20 and not 2020
Corona_Cases.long$Date<-as.Date(Corona_Cases.long$Date,format = "%m.%d.%y")
Corona_Cases.long$Date.numeric<-as.numeric(Corona_Cases.long$Date)
kable(table(select(Corona_Cases.long,c("Country.Region","case_type"))),caption = "Number of death and total case longitudinal datapoints per geographical region")| death | total | |
|---|---|---|
| Afghanistan | 115 | 115 |
| Albania | 115 | 115 |
| Algeria | 115 | 115 |
| Andorra | 115 | 115 |
| Angola | 115 | 115 |
| Antigua and Barbuda | 115 | 115 |
| Argentina | 115 | 115 |
| Armenia | 115 | 115 |
| Australia | 920 | 920 |
| Austria | 115 | 115 |
| Azerbaijan | 115 | 115 |
| Bahamas | 115 | 115 |
| Bahrain | 115 | 115 |
| Bangladesh | 115 | 115 |
| Barbados | 115 | 115 |
| Belarus | 115 | 115 |
| Belgium | 115 | 115 |
| Belize | 115 | 115 |
| Benin | 115 | 115 |
| Bhutan | 115 | 115 |
| Bolivia | 115 | 115 |
| Bosnia and Herzegovina | 115 | 115 |
| Botswana | 115 | 115 |
| Brazil | 115 | 115 |
| Brunei | 115 | 115 |
| Bulgaria | 115 | 115 |
| Burkina Faso | 115 | 115 |
| Burma | 115 | 115 |
| Burundi | 115 | 115 |
| Cabo Verde | 115 | 115 |
| Cambodia | 115 | 115 |
| Cameroon | 115 | 115 |
| Canada | 1610 | 1610 |
| Central African Republic | 115 | 115 |
| Chad | 115 | 115 |
| Chile | 115 | 115 |
| China | 3795 | 3795 |
| Colombia | 115 | 115 |
| Comoros | 115 | 115 |
| Congo (Brazzaville) | 115 | 115 |
| Congo (Kinshasa) | 115 | 115 |
| Costa Rica | 115 | 115 |
| Cote d’Ivoire | 115 | 115 |
| Croatia | 115 | 115 |
| Cuba | 115 | 115 |
| Cyprus | 115 | 115 |
| Czechia | 115 | 115 |
| Denmark | 345 | 345 |
| Diamond Princess | 115 | 115 |
| Djibouti | 115 | 115 |
| Dominica | 115 | 115 |
| Dominican Republic | 115 | 115 |
| Ecuador | 115 | 115 |
| Egypt | 115 | 115 |
| El Salvador | 115 | 115 |
| Equatorial Guinea | 115 | 115 |
| Eritrea | 115 | 115 |
| Estonia | 115 | 115 |
| Eswatini | 115 | 115 |
| Ethiopia | 115 | 115 |
| Fiji | 115 | 115 |
| Finland | 115 | 115 |
| France | 1265 | 1265 |
| Gabon | 115 | 115 |
| Gambia | 115 | 115 |
| Georgia | 115 | 115 |
| Germany | 115 | 115 |
| Ghana | 115 | 115 |
| Greece | 115 | 115 |
| Grenada | 115 | 115 |
| Guatemala | 115 | 115 |
| Guinea | 115 | 115 |
| Guinea-Bissau | 115 | 115 |
| Guyana | 115 | 115 |
| Haiti | 115 | 115 |
| Holy See | 115 | 115 |
| Honduras | 115 | 115 |
| Hungary | 115 | 115 |
| Iceland | 115 | 115 |
| India | 115 | 115 |
| Indonesia | 115 | 115 |
| Iran | 115 | 115 |
| Iraq | 115 | 115 |
| Ireland | 115 | 115 |
| Israel | 115 | 115 |
| Italy | 115 | 115 |
| Jamaica | 115 | 115 |
| Japan | 115 | 115 |
| Jordan | 115 | 115 |
| Kazakhstan | 115 | 115 |
| Kenya | 115 | 115 |
| Korea, South | 115 | 115 |
| Kosovo | 115 | 115 |
| Kuwait | 115 | 115 |
| Kyrgyzstan | 115 | 115 |
| Laos | 115 | 115 |
| Latvia | 115 | 115 |
| Lebanon | 115 | 115 |
| Lesotho | 115 | 115 |
| Liberia | 115 | 115 |
| Libya | 115 | 115 |
| Liechtenstein | 115 | 115 |
| Lithuania | 115 | 115 |
| Luxembourg | 115 | 115 |
| Madagascar | 115 | 115 |
| Malawi | 115 | 115 |
| Malaysia | 115 | 115 |
| Maldives | 115 | 115 |
| Mali | 115 | 115 |
| Malta | 115 | 115 |
| Mauritania | 115 | 115 |
| Mauritius | 115 | 115 |
| Mexico | 115 | 115 |
| Moldova | 115 | 115 |
| Monaco | 115 | 115 |
| Mongolia | 115 | 115 |
| Montenegro | 115 | 115 |
| Morocco | 115 | 115 |
| Mozambique | 115 | 115 |
| MS Zaandam | 115 | 115 |
| Namibia | 115 | 115 |
| Nepal | 115 | 115 |
| Netherlands | 575 | 575 |
| New Zealand | 115 | 115 |
| Nicaragua | 115 | 115 |
| Niger | 115 | 115 |
| Nigeria | 115 | 115 |
| North Macedonia | 115 | 115 |
| Norway | 115 | 115 |
| Oman | 115 | 115 |
| Pakistan | 115 | 115 |
| Panama | 115 | 115 |
| Papua New Guinea | 115 | 115 |
| Paraguay | 115 | 115 |
| Peru | 115 | 115 |
| Philippines | 115 | 115 |
| Poland | 115 | 115 |
| Portugal | 115 | 115 |
| Qatar | 115 | 115 |
| Romania | 115 | 115 |
| Russia | 115 | 115 |
| Rwanda | 115 | 115 |
| Saint Kitts and Nevis | 115 | 115 |
| Saint Lucia | 115 | 115 |
| Saint Vincent and the Grenadines | 115 | 115 |
| San Marino | 115 | 115 |
| Sao Tome and Principe | 115 | 115 |
| Saudi Arabia | 115 | 115 |
| Senegal | 115 | 115 |
| Serbia | 115 | 115 |
| Seychelles | 115 | 115 |
| Sierra Leone | 115 | 115 |
| Singapore | 115 | 115 |
| Slovakia | 115 | 115 |
| Slovenia | 115 | 115 |
| Somalia | 115 | 115 |
| South Africa | 115 | 115 |
| South Sudan | 115 | 115 |
| Spain | 115 | 115 |
| Sri Lanka | 115 | 115 |
| Sudan | 115 | 115 |
| Suriname | 115 | 115 |
| Sweden | 115 | 115 |
| Switzerland | 115 | 115 |
| Syria | 115 | 115 |
| Taiwan* | 115 | 115 |
| Tajikistan | 115 | 115 |
| Tanzania | 115 | 115 |
| Thailand | 115 | 115 |
| Timor-Leste | 115 | 115 |
| Togo | 115 | 115 |
| Trinidad and Tobago | 115 | 115 |
| Tunisia | 115 | 115 |
| Turkey | 115 | 115 |
| Uganda | 115 | 115 |
| Ukraine | 115 | 115 |
| United Arab Emirates | 115 | 115 |
| United Kingdom | 1265 | 1265 |
| Uruguay | 115 | 115 |
| US | 115 | 115 |
| US_state | 375015 | 375015 |
| Uzbekistan | 115 | 115 |
| Venezuela | 115 | 115 |
| Vietnam | 115 | 115 |
| West Bank and Gaza | 115 | 115 |
| Western Sahara | 115 | 115 |
| Yemen | 115 | 115 |
| Zambia | 115 | 115 |
| Zimbabwe | 115 | 115 |
# Decouple population and lat/long data, refactor to make it more tidy
metadata_columns<-c("Lat","Long","Population")
metadata<-unique(select(filter(Corona_Cases.long,case_type=="death"),c("Country.Region","Province.State","City",all_of(metadata_columns))))
Corona_Cases.long<-select(Corona_Cases.long,-all_of(metadata_columns))
# Some counties are not summarized on the country level. collapse all but US
Corona_Cases.long<-rbind.fill(ddply(filter(Corona_Cases.long,!Country.Region=="US_state"),c("case_type","Country.Region","Date","Date.numeric"),summarise,Total_confirmed_cases=sum(Total_confirmed_cases)),filter(Corona_Cases.long,Country.Region=="US_state"))
# Put total case and deaths side-by-side (wide)
Corona_Cases<-spread(Corona_Cases.long,key = case_type,value = Total_confirmed_cases)
#Compute moratlity rate
Corona_Cases$mortality_rate<-Corona_Cases$death/Corona_Cases$total
#TMP
Corona_Cases<-plyr::rename(Corona_Cases,c("total"="Total_confirmed_cases","death"="Total_confirmed_deaths"))
##------------------------------------------
## log10 transform total # cases
##------------------------------------------
Corona_Cases$Total_confirmed_cases.log<-log(Corona_Cases$Total_confirmed_cases,10)
Corona_Cases$Total_confirmed_deaths.log<-log(Corona_Cases$Total_confirmed_deaths,10)
##------------------------------------------
##------------------------------------------
## Compute # of days since 100th for US data
##------------------------------------------
# Find day that 100th case was found for Country/Province. NOTE: Non US countries may have weird provinces. For example, Frane is summairzed at the country level but also had 3 providences. I've only ensured the U.S. case100 works... so the case100_date for U.S. is summarized both for the entire country (regardless of state) and on a per-state level.
# TODO: consider city-level summary as well. This data may be sparse
Corona_Cases<-merge(Corona_Cases,ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,case100_date=min(Date.numeric)))
Corona_Cases$Days_since_100<-Corona_Cases$Date.numeric-Corona_Cases$case100_date
##------------------------------------------
## Add population and lat/long data (CURRENTLY US ONLY)
##------------------------------------------
kable(filter(metadata,(is.na(Country.Region) | is.na(Population) )) %>% select(c("Country.Region","Province.State","City")) %>% unique(),caption = "Regions for which either population or Country is NA")| Country.Region | Province.State | City |
|---|
# Drop missing data
metadata<-filter(metadata,!(is.na(Country.Region) | is.na(Population) ))
# Convert remaining pop to numeric
metadata$Population<-as.numeric(metadata$Population)## Warning: NAs introduced by coercion# Add metadata to cases
Corona_Cases<-merge(Corona_Cases,metadata,all.x = T)
##------------------------------------------
## Compute total and death cases relative to population
##------------------------------------------
Corona_Cases$Total_confirmed_cases.per100<-100*Corona_Cases$Total_confirmed_cases/Corona_Cases$Population
Corona_Cases$Total_confirmed_deaths.per100<-100*Corona_Cases$Total_confirmed_deaths/Corona_Cases$Population
##------------------------------------------
## Filter df for US state-wide stats
##------------------------------------------
Corona_Cases.US_state<-filter(Corona_Cases,Country.Region=="US_state" & Total_confirmed_cases>0 )
kable(table(select(Corona_Cases.US_state,c("Province.State"))),caption = "Number of longitudinal datapoints (total/death) per state")| Var1 | Freq |
|---|---|
| Alabama | 3529 |
| Alaska | 627 |
| Arizona | 911 |
| Arkansas | 3630 |
| California | 3490 |
| Colorado | 3188 |
| Connecticut | 538 |
| Delaware | 219 |
| Diamond Princess | 60 |
| District of Columbia | 61 |
| Florida | 3818 |
| Georgia | 8400 |
| Grand Princess | 61 |
| Guam | 61 |
| Hawaii | 317 |
| Idaho | 1661 |
| Illinois | 4686 |
| Indiana | 4905 |
| Iowa | 4297 |
| Kansas | 3489 |
| Kentucky | 5184 |
| Louisiana | 3565 |
| Maine | 900 |
| Maryland | 1409 |
| Massachusetts | 931 |
| Michigan | 4251 |
| Minnesota | 3915 |
| Mississippi | 4414 |
| Missouri | 4730 |
| Montana | 1490 |
| Nebraska | 2559 |
| Nevada | 653 |
| New Hampshire | 631 |
| New Jersey | 1384 |
| New Mexico | 1430 |
| New York | 3453 |
| North Carolina | 5139 |
| North Dakota | 1640 |
| Northern Mariana Islands | 46 |
| Ohio | 4597 |
| Oklahoma | 3396 |
| Oregon | 1778 |
| Pennsylvania | 3663 |
| Puerto Rico | 61 |
| Rhode Island | 358 |
| South Carolina | 2593 |
| South Dakota | 2117 |
| Tennessee | 4948 |
| Texas | 10069 |
| Utah | 850 |
| Vermont | 840 |
| Virgin Islands | 61 |
| Virginia | 6422 |
| Washington | 2371 |
| West Virginia | 2268 |
| Wisconsin | 3420 |
| Wyoming | 1060 |
Corona_Cases.US_state<-merge(Corona_Cases.US_state,ddply(filter(Corona_Cases.US_state,Total_confirmed_cases>100),c("Province.State"),summarise,case100_date_state=min(Date.numeric)))
Corona_Cases.US_state$Days_since_100_state<-Corona_Cases.US_state$Date.numeric-Corona_Cases.US_state$case100_date_stateANALYSIS
Q1: What is the trend in cases, mortality across geopgraphical regions?
Plot # of cases vs time
* For each geographical set:
* comparative longitudinal case trend (absolute & log scale)
* comparative longitudinal mortality trend
* death vs total correlation
| question | dataset | x | y | color | facet | pch | dimentions |
|---|---|---|---|---|---|---|---|
| comparative_longitudinal_case_trend | long | time | log_cases | geography | none (case type?) | case_type | [15, 50, 4] geography x (2 scale?) case type |
| comparative longitudinal case trend | long | time | cases | geography | case_type | ? | [15, 50, 4] geography x (2+ scale) case type |
| comparative longitudinal mortality trend | wide | time | mortality rate | geography | none | none | [15, 50, 4] geography |
| death vs total correlation | wide | cases | deaths | geography | none | none | [15, 50, 4] geography |
# total cases vs time
# death cases vs time
# mortality rate vs time
# death vs mortality
# death vs mortality
# total & death case vs time (same plot)
#<question> <x> <y> <colored> <facet> <dataset>
## trend in case/deaths over time, comapred across regions <time> <log cases> <geography*> <none> <.wide>
## trend in case/deaths over time, comapred across regions <time> <cases> <geography*> <case_type> <.long>
## trend in mortality rate over time, comapred across regions <time> <mortality rate> <geography*> <none>
## how are death/mortality related/correlated? <time> <log cases> <geography*> <none>
## how are death and case load correlated? <cases> <deaths>
# lm for each?? - > apply lm from each region starting from 100th case. m, b associated with each.
# input: geographical regsion, logcase vs day (100th case)
# output: m, b for each geographical region ID
#total/death on same plot- diffeer by 2 logs, so when plotting log, use pch. when plotting absolute, need to use free scales
#when plotting death and case on same, melt.
#CoronaCases - > filter sets (3)
#world - choose countries with sufficent data
N<-ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,n=length(Country.Region))
ggplot(filter(N,n<100),aes(x=n))+
geom_histogram()+
default_theme+
ggtitle("Distribution of number of days with at least 100 confirmed cases for each region")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
kable(arrange(N,-n),caption="Sorted number of days with at least 100 confirmed cases")| Country.Region | n |
|---|---|
| US_state | 27987 |
| China | 115 |
| Diamond Princess | 96 |
| Korea, South | 86 |
| Japan | 85 |
| Italy | 83 |
| Iran | 80 |
| Singapore | 77 |
| France | 76 |
| Germany | 76 |
| Spain | 75 |
| US | 74 |
| Switzerland | 72 |
| United Kingdom | 72 |
| Belgium | 71 |
| Netherlands | 71 |
| Norway | 71 |
| Sweden | 71 |
| Austria | 69 |
| Malaysia | 68 |
| Australia | 67 |
| Bahrain | 67 |
| Denmark | 67 |
| Canada | 66 |
| Qatar | 66 |
| Iceland | 65 |
| Brazil | 64 |
| Czechia | 64 |
| Finland | 64 |
| Greece | 64 |
| Iraq | 64 |
| Israel | 64 |
| Portugal | 64 |
| Slovenia | 64 |
| Egypt | 63 |
| Estonia | 63 |
| India | 63 |
| Ireland | 63 |
| Kuwait | 63 |
| Philippines | 63 |
| Poland | 63 |
| Romania | 63 |
| Saudi Arabia | 63 |
| Indonesia | 62 |
| Lebanon | 62 |
| San Marino | 62 |
| Thailand | 62 |
| Chile | 61 |
| Pakistan | 61 |
| Luxembourg | 60 |
| Peru | 60 |
| Russia | 60 |
| Ecuador | 59 |
| Mexico | 59 |
| Slovakia | 59 |
| South Africa | 59 |
| United Arab Emirates | 59 |
| Armenia | 58 |
| Colombia | 58 |
| Croatia | 58 |
| Panama | 58 |
| Serbia | 58 |
| Taiwan* | 58 |
| Turkey | 58 |
| Argentina | 57 |
| Bulgaria | 57 |
| Latvia | 57 |
| Uruguay | 57 |
| Algeria | 56 |
| Costa Rica | 56 |
| Dominican Republic | 56 |
| Hungary | 56 |
| Andorra | 55 |
| Bosnia and Herzegovina | 55 |
| Jordan | 55 |
| Lithuania | 55 |
| Morocco | 55 |
| New Zealand | 55 |
| North Macedonia | 55 |
| Vietnam | 55 |
| Albania | 54 |
| Cyprus | 54 |
| Malta | 54 |
| Moldova | 54 |
| Brunei | 53 |
| Burkina Faso | 53 |
| Sri Lanka | 53 |
| Tunisia | 53 |
| Ukraine | 52 |
| Azerbaijan | 51 |
| Ghana | 51 |
| Kazakhstan | 51 |
| Oman | 51 |
| Senegal | 51 |
| Venezuela | 51 |
| Afghanistan | 50 |
| Cote d’Ivoire | 50 |
| Cuba | 49 |
| Mauritius | 49 |
| Uzbekistan | 49 |
| Cambodia | 48 |
| Cameroon | 48 |
| Honduras | 48 |
| Nigeria | 48 |
| West Bank and Gaza | 48 |
| Belarus | 47 |
| Georgia | 47 |
| Bolivia | 46 |
| Kosovo | 46 |
| Kyrgyzstan | 46 |
| Montenegro | 46 |
| Congo (Kinshasa) | 45 |
| Kenya | 44 |
| Niger | 43 |
| Guinea | 42 |
| Rwanda | 42 |
| Trinidad and Tobago | 42 |
| Paraguay | 41 |
| Bangladesh | 40 |
| Djibouti | 38 |
| El Salvador | 37 |
| Guatemala | 36 |
| Madagascar | 35 |
| Mali | 34 |
| Congo (Brazzaville) | 31 |
| Jamaica | 31 |
| Gabon | 29 |
| Somalia | 29 |
| Tanzania | 29 |
| Ethiopia | 28 |
| Burma | 27 |
| Sudan | 26 |
| Liberia | 25 |
| Maldives | 23 |
| Equatorial Guinea | 22 |
| Cabo Verde | 20 |
| Sierra Leone | 18 |
| Guinea-Bissau | 17 |
| Togo | 17 |
| Zambia | 16 |
| Eswatini | 15 |
| Chad | 14 |
| Tajikistan | 13 |
| Haiti | 11 |
| Sao Tome and Principe | 11 |
| Benin | 9 |
| Nepal | 9 |
| Uganda | 9 |
| Central African Republic | 8 |
| South Sudan | 8 |
| Guyana | 6 |
| Mozambique | 5 |
| Yemen | 1 |
# Pick top 15 countries with data
max_colors<-12
# find way to fix this- China has diff provences. Plot doesnt look right...
sufficient_data<-arrange(filter(N,!Country.Region %in% c("US_state", "Diamond Princess")),-n)[1:max_colors,]
kable(sufficient_data,caption = paste0("Top ",max_colors," countries with sufficient data"))| Country.Region | n |
|---|---|
| China | 115 |
| Korea, South | 86 |
| Japan | 85 |
| Italy | 83 |
| Iran | 80 |
| Singapore | 77 |
| France | 76 |
| Germany | 76 |
| Spain | 75 |
| US | 74 |
| Switzerland | 72 |
| United Kingdom | 72 |
Corona_Cases.world<-filter(Corona_Cases,Country.Region %in% c(sufficient_data$Country.Region))
#us
# - by state
Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
# summarize
#!City %in% c("Unassigned")
# - specific cities
#mortality_rate!=Inf & mortality_rate<=1
Corona_Cases.UScity<-filter(Corona_Cases,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") & City %in% c("Bucks","Baltimore City", "New York","Burlington"))
measure_vars_long<-c("Total_confirmed_cases.log","Total_confirmed_cases","Total_confirmed_deaths","Total_confirmed_deaths.log")
melt_arg_list<-list(variable.name = "case_type",value.name = "cases",measure.vars = c("Total_confirmed_cases","Total_confirmed_deaths"))
melt_arg_list$data=NULL
melt_arg_list$data=select(Corona_Cases.world,-ends_with(match = "log"))
Corona_Cases.world.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.UScity,-ends_with(match = "log"))
Corona_Cases.UScity.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.US_state,-ends_with(match = "log"))
Corona_Cases.US_state.long<-do.call(melt,melt_arg_list)
Corona_Cases.world.long$cases.log<-log(Corona_Cases.world.long$cases,10)
Corona_Cases.US_state.long$cases.log<-log(Corona_Cases.US_state.long$cases,10)
Corona_Cases.UScity.long$cases.log<-log(Corona_Cases.UScity.long$cases,10)
# what is the current death and total case load for US? For world? For states?
#-absolute
#-log
# what is mortality rate (US, world)
#-absolute
#how is death and case correlated? (US, world)
#-absolute#Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
#Corona_Cases.US.case100<-filter(Corona_Cases.US, Days_since_100>=0)
# linear model parameters
#(model_fit<-lm(formula = Total_confirmed_cases.log~Days_since_100,data= Corona_Cases.US.case100 ))
#(slope<-model_fit$coefficients[2])
#(intercept<-model_fit$coefficients[1])
# Correlation coefficient
#cor(x = Corona_Cases.US.case100$Days_since_100,y = Corona_Cases.US.case100$Total_confirmed_cases.log)
##------------------------------------------
## Plot World Data
##------------------------------------------
# Timestamp for world
timestamp_plot.world<-paste("Most recent date for which data available:",max(Corona_Cases.world$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")
# Base template for plots
baseplot.world<-ggplot(data=NULL,aes(x=Days_since_100,col=Country.Region))+
default_theme+
scale_color_brewer(type = "qualitative",palette = "Paired")+
ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
theme(legend.position = "bottom",plot.title = element_text(size=12))
##/////////////////////////
### Plot Longitudinal cases
(Corona_Cases.world.long.plot<-baseplot.world+
geom_point(data=Corona_Cases.world.long,aes(y=cases))+
geom_line(data=Corona_Cases.world.long,aes(y=cases))+
facet_wrap(~case_type,scales = "free_y",ncol=1)+
ggtitle(timestamp_plot.world)
)
(Corona_Cases.world.loglong.plot<-baseplot.world+
geom_point(data=Corona_Cases.world.long,aes(y=cases.log))+
geom_line(data=Corona_Cases.world.long,aes(y=cases.log))+
facet_wrap(~case_type,scales = "free_y",ncol=1)+
ggtitle(timestamp_plot.world))
##/////////////////////////
### Plot Longitudinal mortality rate
(Corona_Cases.world.mortality.plot<-baseplot.world+
geom_point(data=Corona_Cases.world,aes(y=mortality_rate))+
geom_line(data=Corona_Cases.world,aes(y=mortality_rate))+
ylim(c(0,0.3))+
ggtitle(timestamp_plot.world))## Warning: Removed 100 rows containing missing values (geom_point).## Warning: Removed 100 row(s) containing missing values (geom_path).
##/////////////////////////
### Plot death vs total case correlation
(Corona_Cases.world.casecor.plot<-ggplot(Corona_Cases.world,aes(x=Total_confirmed_cases,y=Total_confirmed_deaths,col=Country.Region))+
geom_point()+
geom_line()+
default_theme+
scale_color_brewer(type = "qualitative",palette = "Paired")+
ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
theme(legend.position = "bottom",plot.title = element_text(size=12))+
ggtitle(timestamp_plot.world))
### Write polots
write_plot(Corona_Cases.world.long.plot,wd = results_dir)## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.long.plot.png"write_plot(Corona_Cases.world.loglong.plot,wd = results_dir)## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.loglong.plot.png"write_plot(Corona_Cases.world.mortality.plot,wd = results_dir)## Warning: Removed 100 rows containing missing values (geom_point).
## Warning: Removed 100 row(s) containing missing values (geom_path).## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.mortality.plot.png"write_plot(Corona_Cases.world.casecor.plot,wd = results_dir)## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.casecor.plot.png"##------------------------------------------
## Plot US State Data
##-----------------------------------------
baseplot.US<-ggplot(data=NULL,aes(x=Days_since_100_state,col=case_type))+
default_theme+
facet_wrap(~Province.State)+
ggtitle(paste("Log10 cases over time,",timestamp_plot.world))
Corona_Cases.US_state.long.plot<-baseplot.US+geom_point(data=Corona_Cases.US_state.long,aes(y=cases.log))
##------------------------------------------
## Plot US City Data
##-----------------------------------------
Corona_Cases.US.plotdata<-filter(Corona_Cases.US_state,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") &
City %in% c("Bucks","Baltimore City", "New York","Burlington") &
Total_confirmed_cases>0)
timestamp_plot<-paste("Most recent date for which data available:",max(Corona_Cases.US.plotdata$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")
city_colors<-c("Bucks"='#beaed4',"Baltimore City"='#386cb0', "New York"='#7fc97f',"Burlington"='#fdc086')
##/////////////////////////
### Plot death vs total case correlation
(Corona_Cases.city.loglong.plot<-ggplot(melt(Corona_Cases.US.plotdata,measure.vars = c("Total_confirmed_cases.log","Total_confirmed_deaths.log"),variable.name = "case_type",value.name = "cases"),aes(x=Date,y=cases,col=City,pch=case_type))+
geom_point(size=4)+
geom_line()+
default_theme+
#facet_wrap(~case_type)+
ggtitle(paste("Log10 total and death cases over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
scale_color_manual(values = city_colors)+
scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))
(Corona_Cases.city.long.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases,col=City))+
geom_point(size=4)+
geom_line()+
default_theme+
facet_grid(~Province.State,scales = "free_y")+
ggtitle(paste("MD, PA, NJ total cases over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))
+
scale_color_manual(values = city_colors)+
scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))
(Corona_Cases.city.mortality.plot<-ggplot(Corona_Cases.US.plotdata,aes(x=Date,y=mortality_rate,col=City))+
geom_point(size=3)+
geom_line(size=2)+
default_theme+
ggtitle(paste("Mortality rate (deaths/total) over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
scale_color_manual(values = city_colors)+
scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))
(Corona_Cases.city.casecor.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(y=Total_confirmed_deaths,x=Total_confirmed_cases,col=City))+
geom_point(size=3)+
geom_line(size=2)+
default_theme+
ggtitle(paste("Correlation of death vs total cases,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12))+
scale_color_manual(values = city_colors))
(Corona_Cases.city.long.normalized.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases.per100,col=City))+
geom_point(size=4)+
geom_line()+
default_theme+
facet_grid(~Province.State)+
ggtitle(paste("MD, PA, NJ total cases over time per 100 people,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
scale_color_manual(values = city_colors) +
scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))
write_plot(Corona_Cases.city.long.plot,wd = results_dir_custom)## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.plot.png"write_plot(Corona_Cases.city.loglong.plot,wd = results_dir_custom)## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.loglong.plot.png"write_plot(Corona_Cases.city.mortality.plot,wd = results_dir_custom)## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.mortality.plot.png"write_plot(Corona_Cases.city.casecor.plot,wd = results_dir_custom)## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.casecor.plot.png"write_plot(Corona_Cases.city.long.normalized.plot,wd = results_dir_custom)## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.normalized.plot.png"Q1b what is the model
Fit the cases to a linear model 1. Find time at which the case vs date becomes linear in each plot
2. Fit linear model for each city
# What is the predict # of cases for the next few days?
# How is the model performing historically?
Corona_Cases.US_state.summary<-ddply(Corona_Cases.US_state,
c("Province.State","Date"),
summarise,
Total_confirmed_cases_perstate=sum(Total_confirmed_cases)) %>%
filter(Total_confirmed_cases_perstate>100)
# Compute the states with the most cases (for coloring and for linear model)
top_states_totals<-head(ddply(Corona_Cases.US_state.summary,c("Province.State"),summarise, Total_confirmed_cases_perstate.max=max(Total_confirmed_cases_perstate)) %>% arrange(-Total_confirmed_cases_perstate.max),n=max_colors)
kable(top_states_totals,caption = "Top 12 States, total count ")
top_states<-top_states_totals$Province.State
# Manually fix states so that Maryland is switched out for New York
top_states_modified<-c(top_states[top_states !="New York"],"Maryland")
# Plot with all states:
(Corona_Cases.US_state.summary.plot<-ggplot(Corona_Cases.US_state.summary,aes(x=Date,y=Total_confirmed_cases_perstate))+
geom_point()+
geom_point(data=filter(Corona_Cases.US_state.summary,Province.State %in% top_states),aes(col=Province.State))+
scale_color_brewer(type = "qualitative",palette = "Paired")+
default_theme+
theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
ggtitle("Total confirmed cases per state, top 12 colored")+
scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))
##------------------------------------------
## Fit linear model to time vs total cases
##-----------------------------------------
# First, find the date at which each state's cases vs time becomes lienar (2nd derivative is about 0)
li<-ddply(Corona_Cases.US_state.summary,c("Province.State"),find_linear_index)
# Compute linear model for each state starting at the point at which data becomes linear
for(i in 1:nrow(li)){
Province.State.i<-li[i,"Province.State"]
date.i<-li[i,"V1"]
data.i<-filter(Corona_Cases.US_state.summary,Province.State==Province.State.i & as.numeric(Date) >= date.i)
model_results<-lm(data.i,formula = Total_confirmed_cases_perstate~Date)
slope<-model_results$coefficients[2]
intercept<-model_results$coefficients[1]
li[li$Province.State==Province.State.i,"m"]<-slope
li[li$Province.State==Province.State.i,"b"]<-intercept
}
# Compute top state case load with fitted model
(Corona_Cases.US_state.lm.plot<-ggplot(filter(Corona_Cases.US_state.summary,Province.State %in% top_states_modified ))+
geom_abline(data=filter(li,Province.State %in% top_states_modified),
aes(slope = m,intercept = b,col=Province.State),lty=2)+
geom_point(aes(x=Date,y=Total_confirmed_cases_perstate,col=Province.State))+
scale_color_brewer(type = "qualitative",palette = "Paired")+
default_theme+
theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
ggtitle("Total confirmed cases per state, top 12 colored")+
scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))
##------------------------------------------
## Predict the number of total cases over the next week
##-----------------------------------------
predicted_days<-c(0,1,2,3,7)+as.numeric(as.Date("2020-04-20"))
predicted_days_df<-data.frame(matrix(ncol=3))
names(predicted_days_df)<-c("Province.State","days","Total_confirmed_cases_perstate")
# USe model parameters to estiamte case loads
for(state.i in top_states_modified){
predicted_days_df<-rbind(predicted_days_df,
data.frame(Province.State=state.i,
prediction_model(m = li[li$Province.State==state.i,"m"],
b =li[li$Province.State==state.i,"b"] ,
days =predicted_days )))
}
predicted_days_df$Date<-as.Date(predicted_days_df$days,origin="1970-01-01")
kable(predicted_days_df,caption = "Predicted total cases over the next week for selected states")
##------------------------------------------
## Write plots
##-----------------------------------------
write_plot(Corona_Cases.US_state.summary.plot,wd = results_dir)
write_plot(Corona_Cases.US_state.lm.plot,wd = results_dir)
##------------------------------------------
## Write tables
##-----------------------------------------
write.csv(predicted_days_df,file = paste0(results_dir,"predicted_total_cases_days.csv"),quote = F,row.names = F)Q2: What is the predicted number of cases?
What is the prediction of COVID-19 based on model thus far? Additional questions:
WHy did it take to day 40 to start a log linear trend? How long will it be till x number of cases? When will the plateu happen? Are any effects noticed with social distancing? Delays
##------------------------------------------
## Prediction and Prediction Accuracy
##------------------------------------------
today_num<-max(Corona_Cases.US$Days_since_100)
predicted_days<-today_num+c(1,2,3,7)
#mods = dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)
#today:
Corona_Cases.US[Corona_Cases.US$Days_since_100==(today_num-1),]
Corona_Cases.US[Corona_Cases.US$Days_since_100==today_num,]
Corona_Cases.US$type<-"Historical"
#prediction_values<-prediction_model(m=slope,b=intercept,days = predicted_days)$Total_confirmed_cases
histoical_model<-data.frame(date=today_num,m=slope,b=intercept)
tmp<-data.frame(state=rep(c("A","B"),each=3),x=c(1,2,3,4,5,6))
tmp$y<-c(tmp[1:3,"x"]+5,tmp[4:6,"x"]*5+1)
ddply(tmp,c("state"))
lm(data =tmp,formula = y~x )
train_lm<-function(input_data,subset_coulmn,formula_input){
case_models <- dlply(input_data, subset_coulmn, lm, formula = formula_input)
case_models.parameters <- ldply(case_models, coef)
case_models.parameters<-rename(case_models.parameters,c("b"="(Intercept)","m"=subset_coulmn))
return(case_models.parameters)
}
train_lm(tmp,"state")
dlply(input_data, subset_coulmn, lm,m=)
# model for previous y days
#historical_model_predictions<-data.frame(day_x=NULL,Days_since_100=NULL,Total_confirmed_cases=NULL,Total_confirmed_cases.log=NULL)
# for(i in c(1,2,3,4,5,6,7,8,9,10)){
# #i<-1
# day_x<-today_num-i # 1, 2, 3, 4
# day_x_nextweek<-day_x+c(1,2,3)
# model_fit_x<-lm(data = filter(Corona_Cases.US.case100,Days_since_100 < day_x),formula = Total_confirmed_cases.log~Days_since_100)
# prediction_day_x_nextweek<-prediction_model(m = model_fit_x$coefficients[2],b = model_fit_x$coefficients[1],days = day_x_nextweek)
# prediction_day_x_nextweek$type<-"Predicted"
# acutal_day_x_nextweek<-filter(Corona_Cases.US,Days_since_100 %in% day_x_nextweek) %>% select(c(Days_since_100,Total_confirmed_cases,Total_confirmed_cases.log))
# acutal_day_x_nextweek$type<-"Historical"
# historical_model_predictions.i<-data.frame(day_x=day_x,rbind(acutal_day_x_nextweek,prediction_day_x_nextweek))
# historical_model_predictions<-rbind(historical_model_predictions.i,historical_model_predictions)
# }
#historical_model_predictions.withHx<-rbind.fill(historical_model_predictions,data.frame(Corona_Cases.US,type="Historical"))
#historical_model_predictions.withHx$Total_confirmed_cases.log2<-log(historical_model_predictions.withHx$Total_confirmed_cases,2)
(historical_model_predictions.plot<-ggplot(historical_model_predictions.withHx,aes(x=Days_since_100,y=Total_confirmed_cases.log,col=type))+
geom_point(size=3)+
default_theme+
theme(legend.position = "bottom")+
#geom_abline(slope = slope,intercept =intercept,lty=2)+
#facet_wrap(~case_type,ncol=1)+
scale_color_manual(values = c("Historical"="#377eb8","Predicted"="#e41a1c")))
write_plot(historical_model_predictions.plot,wd=results_dir)Q3: What is the effect on social distancing, descreased mobility on case load?
Load data from Google which compoutes % change in user mobility relative to baseline for * Recreation
* Workplace
* Residence
* Park
* Grocery
Data from https://www.google.com/covid19/mobility/
# See pre-processing section for script on gathering mobility data
# UNDER DEVELOPMENT
mobility<-read.csv("/Users/stevensmith/Projects/MIT_COVID19/mobility.csv",header = T,stringsAsFactors = F)
#mobility$Retail_Recreation<-as.numeric(sub(mobility$Retail_Recreation,pattern = "%",replacement = ""))
#mobility$Workplace<-as.numeric(sub(mobility$Workplace,pattern = "%",replacement = ""))
#mobility$Residential<-as.numeric(sub(mobility$Residential,pattern = "%",replacement = ""))
##------------------------------------------
## Show relationship between mobility and caseload
##------------------------------------------
mobility$County<-gsub(mobility$County,pattern = " County",replacement = "")
Corona_Cases.US_state.mobility<-merge(Corona_Cases.US_state,plyr::rename(mobility,c("State"="Province.State","County"="City")))
#Corona_Cases.US_state.tmp<-merge(metadata,Corona_Cases.US_state.tmp)
# Needs to happen upsteam, see todos
#Corona_Cases.US_state.tmp$Total_confirmed_cases.perperson<-Corona_Cases.US_state.tmp$Total_confirmed_cases/as.numeric(Corona_Cases.US_state.tmp$Population)
mobility_measures<-c("Retail_Recreation","Grocery_Pharmacy","Parks","Transit","Workplace","Residential")
plot_data<-filter(Corona_Cases.US_state.mobility, Date.numeric==max(Corona_Cases.US_state$Date.numeric) ) %>% melt(measure.vars=mobility_measures)
plot_data$value<-as.numeric(gsub(plot_data$value,pattern = "%",replacement = ""))
plot_data<-filter(plot_data,!is.na(value))
(mobility.plot<-ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
facet_grid(Province.State~variable,scales = "free")+
xlab("Mobility change from baseline (%)")+
ylab(paste0("Confirmed cases per 100 people(Today)"))+
default_theme+
ggtitle("Mobility change vs cases"))
(mobility.global.plot<-ggplot(plot_data,aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
facet_wrap(~variable,scales = "free")+
xlab("Mobility change from baseline (%)")+
ylab(paste0("Confirmed cases (Today) per 100 people"))+
default_theme+
ggtitle("Mobility change vs cases"))
plot_data.permobility_summary<-ddply(plot_data,c("Province.State","variable"),summarise,cor=cor(y =Total_confirmed_cases.per100,x=value),median_change=median(x=value)) %>% arrange(-abs(cor))
kable(plot_data.permobility_summary,caption = "Ranked per-state mobility correlation with total confirmed cases")| Province.State | variable | cor | median_change |
|---|---|---|---|
| Alaska | Transit | -1.0000000 | -63.0 |
| Delaware | Retail_Recreation | 1.0000000 | -39.5 |
| Delaware | Grocery_Pharmacy | 1.0000000 | -17.5 |
| Delaware | Parks | -1.0000000 | 20.5 |
| Delaware | Transit | 1.0000000 | -37.0 |
| Delaware | Workplace | 1.0000000 | -37.0 |
| Delaware | Residential | -1.0000000 | 14.0 |
| Hawaii | Grocery_Pharmacy | 0.9994339 | -34.0 |
| Hawaii | Retail_Recreation | 0.9952917 | -56.0 |
| New Hampshire | Parks | 0.9538409 | -20.0 |
| Connecticut | Grocery_Pharmacy | -0.9044878 | -6.0 |
| Alaska | Residential | 0.8967338 | 13.0 |
| South Dakota | Parks | 0.8923021 | -26.0 |
| Maine | Transit | -0.8904520 | -50.0 |
| Vermont | Parks | 0.8570246 | -35.5 |
| Utah | Workplace | -0.8514789 | -37.0 |
| Wyoming | Transit | -0.8493190 | -16.0 |
| Hawaii | Parks | 0.8211293 | -72.0 |
| Utah | Residential | -0.8169973 | 12.0 |
| Alaska | Grocery_Pharmacy | -0.8134720 | -7.0 |
| Utah | Parks | -0.7907146 | 17.0 |
| Connecticut | Transit | -0.7899244 | -50.0 |
| Massachusetts | Workplace | -0.7866514 | -39.0 |
| Hawaii | Transit | 0.7714358 | -89.0 |
| Utah | Transit | -0.7582161 | -18.0 |
| Rhode Island | Workplace | -0.7474680 | -39.5 |
| Alaska | Workplace | -0.7214938 | -34.0 |
| Maine | Workplace | -0.6746828 | -30.0 |
| Vermont | Grocery_Pharmacy | -0.6516949 | -25.0 |
| New Jersey | Workplace | -0.6473653 | -44.0 |
| New York | Workplace | -0.6414487 | -34.5 |
| Hawaii | Residential | -0.6363072 | 19.0 |
| Montana | Workplace | -0.6351999 | -40.5 |
| Rhode Island | Residential | -0.6263357 | 18.5 |
| Rhode Island | Retail_Recreation | -0.6088746 | -45.0 |
| Arizona | Grocery_Pharmacy | -0.5927171 | -15.0 |
| Nebraska | Workplace | 0.5894447 | -32.0 |
| New York | Retail_Recreation | -0.5836330 | -46.0 |
| North Dakota | Retail_Recreation | -0.5802968 | -43.5 |
| New Jersey | Parks | -0.5745756 | -6.0 |
| Connecticut | Residential | 0.5667860 | 14.0 |
| Massachusetts | Retail_Recreation | -0.5516475 | -44.0 |
| New Jersey | Retail_Recreation | -0.5503955 | -62.5 |
| Utah | Retail_Recreation | -0.5407251 | -40.0 |
| West Virginia | Parks | 0.5361148 | -33.0 |
| Maine | Parks | 0.5279796 | -31.0 |
| New York | Parks | 0.5195797 | 20.0 |
| Arkansas | Parks | -0.5070562 | -12.0 |
| Wyoming | Workplace | -0.5033621 | -31.0 |
| Connecticut | Workplace | -0.5009838 | -39.0 |
| Nebraska | Residential | -0.4934459 | 14.0 |
| Montana | Parks | -0.4913929 | -58.0 |
| New Jersey | Grocery_Pharmacy | -0.4857206 | 2.5 |
| Wisconsin | Transit | -0.4857099 | -23.5 |
| Iowa | Parks | -0.4846011 | 28.5 |
| Arizona | Retail_Recreation | -0.4832532 | -42.5 |
| Idaho | Workplace | -0.4815937 | -29.0 |
| Rhode Island | Parks | 0.4782766 | 52.0 |
| Connecticut | Retail_Recreation | -0.4758725 | -45.0 |
| New Mexico | Grocery_Pharmacy | -0.4738477 | -11.0 |
| New Mexico | Parks | 0.4719227 | -31.5 |
| Montana | Residential | 0.4679731 | 14.0 |
| Illinois | Transit | -0.4576885 | -31.0 |
| Massachusetts | Grocery_Pharmacy | -0.4571418 | -7.0 |
| Wyoming | Parks | -0.4546342 | -4.0 |
| Idaho | Transit | -0.4534782 | -30.0 |
| California | Transit | -0.4515787 | -42.0 |
| California | Residential | 0.4477273 | 14.0 |
| New Mexico | Residential | 0.4455109 | 13.5 |
| New Jersey | Transit | -0.4356006 | -50.5 |
| Kentucky | Parks | -0.4334034 | 28.5 |
| Montana | Retail_Recreation | -0.4329731 | -51.0 |
| Pennsylvania | Workplace | -0.4306945 | -36.0 |
| Vermont | Residential | 0.4282852 | 11.5 |
| New Hampshire | Residential | -0.4099784 | 14.0 |
| South Carolina | Workplace | 0.3964847 | -30.0 |
| Idaho | Grocery_Pharmacy | -0.3947512 | -4.5 |
| Montana | Grocery_Pharmacy | -0.3906632 | -16.0 |
| Alabama | Workplace | -0.3897406 | -29.0 |
| Montana | Transit | -0.3892824 | -41.0 |
| Alabama | Transit | -0.3826796 | -36.5 |
| New Mexico | Retail_Recreation | -0.3786975 | -42.0 |
| Maryland | Grocery_Pharmacy | -0.3778669 | -10.0 |
| New York | Transit | -0.3742053 | -48.0 |
| North Dakota | Workplace | 0.3730090 | -38.0 |
| Maryland | Workplace | -0.3714647 | -35.0 |
| Florida | Residential | 0.3704070 | 14.0 |
| California | Parks | -0.3612644 | -38.5 |
| Pennsylvania | Retail_Recreation | -0.3535914 | -45.0 |
| Hawaii | Workplace | 0.3487421 | -46.0 |
| Idaho | Retail_Recreation | -0.3454422 | -40.5 |
| Arizona | Residential | 0.3440607 | 13.0 |
| California | Retail_Recreation | -0.3393992 | -44.0 |
| Nebraska | Grocery_Pharmacy | 0.3391541 | -0.5 |
| Kansas | Parks | 0.3382521 | 72.0 |
| Nevada | Transit | -0.3319415 | -20.0 |
| Alabama | Grocery_Pharmacy | -0.3284394 | -2.0 |
| Arizona | Transit | 0.3253041 | -38.0 |
| California | Workplace | -0.3210955 | -36.0 |
| Vermont | Retail_Recreation | 0.3204559 | -57.0 |
| Wyoming | Grocery_Pharmacy | -0.3190166 | -9.0 |
| Minnesota | Transit | -0.3188445 | -28.5 |
| Pennsylvania | Parks | 0.3159243 | 13.0 |
| Alaska | Retail_Recreation | 0.3155306 | -39.0 |
| California | Grocery_Pharmacy | -0.3141003 | -12.0 |
| Florida | Parks | -0.3120053 | -43.0 |
| Colorado | Residential | 0.3020027 | 14.0 |
| Michigan | Parks | 0.2994274 | 30.0 |
| Texas | Workplace | 0.2965868 | -32.0 |
| Nevada | Retail_Recreation | -0.2914979 | -43.0 |
| Arkansas | Retail_Recreation | -0.2896251 | -30.0 |
| Maine | Retail_Recreation | -0.2885308 | -42.0 |
| Kansas | Workplace | 0.2877426 | -33.0 |
| Texas | Residential | -0.2872139 | 15.0 |
| Rhode Island | Grocery_Pharmacy | 0.2822096 | -7.5 |
| Rhode Island | Transit | -0.2807816 | -56.0 |
| North Carolina | Grocery_Pharmacy | 0.2787648 | 0.0 |
| Utah | Grocery_Pharmacy | -0.2772232 | -4.0 |
| North Dakota | Parks | 0.2768748 | -34.0 |
| Illinois | Workplace | -0.2757077 | -30.0 |
| Tennessee | Residential | 0.2727332 | 11.5 |
| Texas | Parks | 0.2705256 | -42.0 |
| Tennessee | Workplace | -0.2681764 | -31.0 |
| Virginia | Transit | -0.2669523 | -33.0 |
| Maryland | Retail_Recreation | -0.2666827 | -39.0 |
| Oregon | Grocery_Pharmacy | 0.2661213 | -7.0 |
| West Virginia | Grocery_Pharmacy | -0.2592091 | -6.0 |
| Wisconsin | Parks | 0.2557439 | 51.5 |
| Nevada | Residential | 0.2552158 | 17.0 |
| Georgia | Grocery_Pharmacy | -0.2550690 | -10.0 |
| North Dakota | Grocery_Pharmacy | -0.2498231 | -8.0 |
| North Carolina | Transit | 0.2491852 | -32.0 |
| Michigan | Workplace | -0.2454241 | -40.0 |
| Maryland | Residential | 0.2424859 | 15.0 |
| North Carolina | Workplace | 0.2422922 | -31.5 |
| Alabama | Parks | 0.2409578 | -1.0 |
| New York | Grocery_Pharmacy | -0.2344692 | 8.0 |
| Mississippi | Residential | 0.2330838 | 13.0 |
| North Carolina | Residential | 0.2313008 | 13.0 |
| Pennsylvania | Grocery_Pharmacy | -0.2310886 | -6.0 |
| West Virginia | Workplace | 0.2255598 | -33.0 |
| Iowa | Transit | 0.2241636 | -24.0 |
| South Dakota | Transit | -0.2227811 | -40.0 |
| Maine | Grocery_Pharmacy | -0.2216457 | -13.0 |
| New Jersey | Residential | 0.2192392 | 18.0 |
| Arkansas | Residential | 0.2182307 | 12.0 |
| Oregon | Residential | 0.2126030 | 10.5 |
| Mississippi | Grocery_Pharmacy | -0.2061887 | -8.0 |
| Washington | Workplace | -0.2046620 | -38.0 |
| Kansas | Grocery_Pharmacy | -0.2038763 | -14.0 |
| Nebraska | Parks | 0.2026480 | 55.5 |
| Georgia | Workplace | -0.1990484 | -33.5 |
| Oklahoma | Parks | -0.1967473 | -18.5 |
| Illinois | Parks | 0.1949071 | 26.5 |
| South Carolina | Parks | -0.1932846 | -23.0 |
| Colorado | Parks | -0.1925603 | 2.0 |
| Virginia | Grocery_Pharmacy | -0.1923888 | -8.0 |
| Texas | Transit | 0.1843264 | -42.0 |
| Georgia | Retail_Recreation | -0.1802440 | -41.0 |
| Washington | Transit | -0.1799578 | -33.5 |
| Vermont | Workplace | -0.1794093 | -43.0 |
| Massachusetts | Transit | -0.1773359 | -45.0 |
| Illinois | Residential | 0.1718200 | 14.0 |
| Wyoming | Retail_Recreation | -0.1711615 | -40.0 |
| Virginia | Parks | 0.1707679 | 6.0 |
| Florida | Workplace | -0.1702159 | -33.0 |
| Maine | Residential | -0.1689406 | 11.0 |
| Michigan | Retail_Recreation | -0.1683489 | -53.0 |
| Ohio | Transit | 0.1679777 | -28.0 |
| New Mexico | Transit | 0.1665808 | -38.0 |
| Oklahoma | Residential | 0.1620659 | 15.0 |
| Virginia | Residential | 0.1592110 | 14.0 |
| Connecticut | Parks | 0.1540410 | 43.0 |
| South Dakota | Retail_Recreation | -0.1515385 | -38.5 |
| Massachusetts | Residential | 0.1507185 | 15.0 |
| Wisconsin | Residential | -0.1496169 | 14.0 |
| West Virginia | Retail_Recreation | 0.1468079 | -38.5 |
| West Virginia | Residential | -0.1454239 | 11.0 |
| Oregon | Parks | 0.1454209 | 16.5 |
| Georgia | Residential | -0.1450633 | 13.0 |
| Idaho | Parks | 0.1434363 | -22.0 |
| Arkansas | Workplace | -0.1433232 | -26.0 |
| Missouri | Residential | -0.1412363 | 13.0 |
| Michigan | Grocery_Pharmacy | -0.1388670 | -11.0 |
| Washington | Residential | 0.1335720 | 13.0 |
| Mississippi | Retail_Recreation | -0.1317420 | -40.0 |
| Mississippi | Workplace | -0.1310166 | -33.0 |
| Indiana | Residential | 0.1308326 | 12.0 |
| South Carolina | Residential | -0.1306945 | 12.0 |
| Indiana | Retail_Recreation | 0.1295402 | -38.0 |
| Pennsylvania | Transit | -0.1273436 | -41.5 |
| New Hampshire | Transit | -0.1258710 | -57.0 |
| Minnesota | Workplace | -0.1257585 | -33.0 |
| New Hampshire | Retail_Recreation | -0.1244880 | -41.0 |
| Maryland | Transit | -0.1234201 | -39.0 |
| Massachusetts | Parks | 0.1233046 | 39.0 |
| Kansas | Transit | -0.1223983 | -26.5 |
| Nebraska | Retail_Recreation | 0.1217207 | -36.0 |
| Kentucky | Transit | 0.1213997 | -31.0 |
| Missouri | Workplace | 0.1209115 | -28.5 |
| Washington | Parks | 0.1183456 | -3.5 |
| Indiana | Parks | -0.1136117 | 29.0 |
| New Hampshire | Grocery_Pharmacy | -0.1135958 | -6.0 |
| Mississippi | Transit | -0.1126123 | -38.5 |
| Ohio | Residential | 0.1125650 | 14.0 |
| Florida | Retail_Recreation | 0.1119654 | -43.0 |
| Arkansas | Transit | 0.1101722 | -27.0 |
| Washington | Grocery_Pharmacy | 0.1085146 | -7.0 |
| Alabama | Retail_Recreation | 0.1054916 | -39.0 |
| Pennsylvania | Residential | 0.1050204 | 15.0 |
| South Dakota | Grocery_Pharmacy | 0.1050079 | -9.0 |
| Wisconsin | Grocery_Pharmacy | 0.1049565 | -1.0 |
| North Carolina | Parks | -0.1029833 | 7.0 |
| Oregon | Retail_Recreation | 0.1014469 | -41.0 |
| Arizona | Workplace | -0.1002348 | -35.0 |
| Kentucky | Grocery_Pharmacy | 0.1000627 | 4.0 |
| Missouri | Transit | -0.0994961 | -24.5 |
| Oregon | Workplace | -0.0963260 | -31.0 |
| Indiana | Workplace | 0.0960233 | -34.0 |
| South Dakota | Residential | 0.0916072 | 15.0 |
| New York | Residential | 0.0903392 | 17.5 |
| Wyoming | Residential | 0.0899662 | 12.5 |
| Iowa | Retail_Recreation | -0.0898198 | -38.0 |
| Tennessee | Transit | -0.0877457 | -32.0 |
| Wisconsin | Workplace | -0.0869245 | -31.0 |
| Virginia | Retail_Recreation | -0.0839419 | -35.0 |
| Georgia | Parks | 0.0817867 | -6.0 |
| Minnesota | Retail_Recreation | 0.0808979 | -40.5 |
| Virginia | Workplace | -0.0803391 | -31.5 |
| Idaho | Residential | -0.0802768 | 11.0 |
| South Carolina | Transit | 0.0776277 | -45.0 |
| Ohio | Retail_Recreation | 0.0755579 | -36.0 |
| Michigan | Residential | 0.0739321 | 15.0 |
| Indiana | Grocery_Pharmacy | -0.0738303 | -5.5 |
| Kentucky | Retail_Recreation | 0.0710620 | -29.0 |
| Iowa | Workplace | -0.0710237 | -30.5 |
| Minnesota | Parks | 0.0709090 | -9.0 |
| Ohio | Parks | -0.0708964 | 67.5 |
| Oklahoma | Grocery_Pharmacy | -0.0703593 | -1.0 |
| Ohio | Grocery_Pharmacy | 0.0698985 | 0.0 |
| Florida | Transit | -0.0679165 | -49.0 |
| Texas | Grocery_Pharmacy | 0.0678460 | -14.0 |
| Colorado | Grocery_Pharmacy | -0.0606591 | -17.0 |
| Colorado | Retail_Recreation | -0.0593733 | -44.0 |
| Colorado | Transit | 0.0591464 | -36.0 |
| Kentucky | Residential | 0.0575774 | 12.0 |
| Oregon | Transit | 0.0543395 | -27.5 |
| New Hampshire | Workplace | 0.0541501 | -37.0 |
| Texas | Retail_Recreation | 0.0505747 | -40.0 |
| North Dakota | Residential | 0.0494212 | 17.0 |
| Michigan | Transit | 0.0488792 | -46.0 |
| Minnesota | Grocery_Pharmacy | 0.0484529 | -6.0 |
| Washington | Retail_Recreation | -0.0478128 | -42.0 |
| Tennessee | Parks | -0.0473195 | 10.5 |
| Illinois | Retail_Recreation | 0.0468368 | -40.0 |
| Kentucky | Workplace | -0.0451521 | -36.0 |
| North Dakota | Transit | 0.0425947 | -48.0 |
| Nevada | Workplace | 0.0424249 | -40.0 |
| Missouri | Retail_Recreation | -0.0422353 | -36.0 |
| South Carolina | Retail_Recreation | -0.0389990 | -35.0 |
| Missouri | Parks | 0.0371589 | 0.0 |
| Minnesota | Residential | 0.0369288 | 17.0 |
| Oklahoma | Retail_Recreation | 0.0368373 | -31.0 |
| West Virginia | Transit | 0.0324533 | -45.0 |
| Illinois | Grocery_Pharmacy | -0.0321335 | 2.0 |
| Missouri | Grocery_Pharmacy | 0.0305992 | 2.0 |
| Alabama | Residential | 0.0288174 | 11.0 |
| Tennessee | Grocery_Pharmacy | 0.0277956 | 6.0 |
| Colorado | Workplace | -0.0265871 | -39.0 |
| Kansas | Residential | -0.0262042 | 13.0 |
| Vermont | Transit | 0.0255787 | -63.0 |
| Georgia | Transit | 0.0255582 | -35.0 |
| Arkansas | Grocery_Pharmacy | -0.0246259 | 3.5 |
| South Carolina | Grocery_Pharmacy | 0.0245347 | 1.0 |
| Iowa | Grocery_Pharmacy | 0.0216077 | 4.0 |
| Nevada | Parks | 0.0193138 | -12.5 |
| Oklahoma | Workplace | 0.0188960 | -31.0 |
| Arizona | Parks | -0.0173979 | -44.5 |
| Kansas | Retail_Recreation | -0.0168793 | -38.0 |
| North Carolina | Retail_Recreation | 0.0161506 | -34.0 |
| South Dakota | Workplace | -0.0149227 | -35.0 |
| Wisconsin | Retail_Recreation | 0.0147475 | -44.5 |
| Ohio | Workplace | -0.0142561 | -35.0 |
| Florida | Grocery_Pharmacy | 0.0124758 | -14.0 |
| Nebraska | Transit | -0.0083056 | -9.0 |
| Mississippi | Parks | 0.0069208 | -25.0 |
| Nevada | Grocery_Pharmacy | 0.0058160 | -12.5 |
| Tennessee | Retail_Recreation | -0.0056511 | -30.0 |
| New Mexico | Workplace | -0.0049185 | -34.0 |
| Indiana | Transit | -0.0036504 | -29.0 |
| Maryland | Parks | -0.0036090 | 27.0 |
| Oklahoma | Transit | 0.0021416 | -26.0 |
| Iowa | Residential | 0.0006223 | 13.0 |
| Alaska | Parks | NA | 29.0 |
| District of Columbia | Retail_Recreation | NA | -69.0 |
| District of Columbia | Grocery_Pharmacy | NA | -28.0 |
| District of Columbia | Parks | NA | -65.0 |
| District of Columbia | Transit | NA | -69.0 |
| District of Columbia | Workplace | NA | -48.0 |
| District of Columbia | Residential | NA | 17.0 |
# sanity check
ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(x=Total_confirmed_cases.per100,fill=variable))+geom_histogram()+
facet_grid(~Province.State)+
default_theme+
theme(legend.position = "bottom")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
write_plot(mobility.plot,wd = results_dir)## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.plot.png"write_plot(mobility.global.plot,wd = results_dir)## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.global.plot.png"(plot_data.permobility_summary.plot<-ggplot(plot_data.permobility_summary,aes(x=variable,y=median_change))+
geom_jitter(size=2,width=.2)+
#geom_jitter(data=plot_data.permobility_summary %>% arrange(-abs(median_change)) %>% head(n=15),aes(col=Province.State),size=2,width=.2)+
default_theme+
ggtitle("Per-Sate Median Change in Mobility")+
xlab("Mobility Meaure")+
ylab("Median Change from Baseline"))
write_plot(plot_data.permobility_summary.plot,wd = results_dir)## [1] "/Users/stevensmith/Projects/coronavirus/results/plot_data.permobility_summary.plot.png"DELIVERABLE MANIFEST
The following link to commited documents pushed to github. These are provided as a convienence, but note this is a manual process. The generation of reports, plots and tables is not coupled to the execution of this markdown. ## Report This report, html & pdf
Plots
github_root<-"https://github.com/sbs87/coronavirus/blob/master/"
plot_handle<-c("Corona_Cases.world.long.plot",
"Corona_Cases.world.loglong.plot",
"Corona_Cases.world.mortality.plot",
"Corona_Cases.world.casecor.plot",
"Corona_Cases.city.long.plot",
"Corona_Cases.city.loglong.plot",
"Corona_Cases.city.mortality.plot",
"Corona_Cases.city.casecor.plot",
"Corona_Cases.city.long.normalized.plot",
"Corona_Cases.US_state.lm.plot",
"Corona_Cases.US_state.summary.plot")
deliverable_manifest<-data.frame(
name=c("World total & death cases, longitudinal",
"World log total & death cases, longitudinal",
"World mortality",
"World total & death cases, correlation",
"City total & death cases, longitudinal",
"City log total & death cases, longitudinal",
"City mortality",
"City total & death cases, correlation",
"City population normalized total & death cases, longitudinal",
"State total cases (select) with linear model, longitudinal",
"State total cases, longitudinal"),
plot_handle=plot_handle,
link=paste0(github_root,"results/",plot_handle,".png")
)
(tmp<-data.frame(row_out=apply(deliverable_manifest,MARGIN = 1,FUN = function(x) paste(x[1],x[2],x[3],sep=" | "))))## row_out
## 1 World total & death cases, longitudinal | Corona_Cases.world.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
## 2 World log total & death cases, longitudinal | Corona_Cases.world.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
## 3 World mortality | Corona_Cases.world.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
## 4 World total & death cases, correlation | Corona_Cases.world.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
## 5 City total & death cases, longitudinal | Corona_Cases.city.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
## 6 City log total & death cases, longitudinal | Corona_Cases.city.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
## 7 City mortality | Corona_Cases.city.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
## 8 City total & death cases, correlation | Corona_Cases.city.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
## 9 City population normalized total & death cases, longitudinal | Corona_Cases.city.long.normalized.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
## 10 State total cases (select) with linear model, longitudinal | Corona_Cases.US_state.lm.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
## 11 State total cases, longitudinal | Corona_Cases.US_state.summary.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.pngrow_out<-apply(tmp, 2, paste, collapse="\t\n")
Tables
CONCLUSION
Overall, the trends of COVID-19 cases is no longer in log-linear phase for world or U.S. (but some regions like MD are still in the log-linear phase). Mortality rate (deaths/confirmed RNA-based cases) is >1%, with a range depending on region. Mobility is not a strong indicator of caseload (U.S. data).
See table below for detailed breakdown.
| Question | Answer |
|---|---|
| What is the effect on social distancing, descreased mobility on case load? |
There is not a strong apparent effect on decreased mobility (work, grocery, retail) or increased mobility (at residence, parks) on number of confirmed cases, either as a country (U.S.) or state level. California appears to have one of the best correlations, but this is a mixed bag |
| What is the trend in cases, mortality across geopgraphical regions? |
The confirmed total casees and mortality is overall log-linear for most countries, with a trailing off beginning for most (inlcuding U.S.). On the state level, NY, NJ, PA starting to trail off; MD is still in log-linear phase. Mortality and case load are highly correlated for NY, NJ, PA, MD. The mortality rate flucutates for a given region, but is about 3% overall. |
Sandbox
# Geographical heatmap!
install.packages("maps")
library(maps)
library
mi_counties <- map_data("county", "pennsylvania") %>%
select(lon = long, lat, group, id = subregion)
head(mi_counties)
ggplot(mi_counties, aes(lon, lat)) +
geom_point(size = .25, show.legend = FALSE) +
coord_quickmap()
mi_counties$cases<-1:2226
name_overlaps(metadata,Corona_Cases.US_state)
tmp<-merge(Corona_Cases.US_state,metadata)
ggplot(filter(tmp,Province.State=="Pennsylvania"), aes(Long, Lat, group = as.factor(City))) +
geom_polygon(aes(fill = Total_confirmed_cases), colour = "grey50") +
coord_quickmap()
ggplot(Corona_Cases.US_state, aes(Long, Lat))+
geom_polygon(aes(fill = Total_confirmed_cases ), color = "white")+
scale_fill_viridis_c(option = "C")
dev.off()
require(maps)
require(viridis)
world_map <- map_data("world")
ggplot(world_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill="lightgray", colour = "white")
head(world_map)
head(Corona_Cases.US_state)
unique(select(world_map,c("region","group"))) %>% filter()
some.eu.countries <- c(
"US"
)
# Retrievethe map data
some.eu.maps <- map_data("world", region = some.eu.countries)
# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region.lab.data <- some.eu.maps %>%
group_by(region) %>%
summarise(long = mean(long), lat = mean(lat))
unique(filter(some.eu.maps,subregion %in% Corona_Cases.US_state$Province.State) %>% select(subregion))
unique(Corona_Cases.US_state$Total_confirmed_cases.log)
ggplot(filter(Corona_Cases.US_state,Date=="2020-04-17") aes(x = Long, y = Lat)) +
geom_polygon(aes( fill = Total_confirmed_cases.log))+
#geom_text(aes(label = region), data = region.lab.data, size = 3, hjust = 0.5)+
#scale_fill_viridis_d()+
#theme_void()+
theme(legend.position = "none")
library("sf")
library("rnaturalearth")
library("rnaturalearthdata")
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
ggplot(data = world) +
geom_sf()
counties <- st_as_sf(map("county", plot = FALSE, fill = TRUE))
counties <- subset(counties, grepl("florida", counties$ID))
counties$area <- as.numeric(st_area(counties))
#install.packages("lwgeom")
class(counties)
head(counties)
ggplot(data = world) +
geom_sf(data=Corona_Cases.US_state) +
#geom_sf(data = counties, aes(fill = area)) +
geom_sf(data = counties, aes(fill = area)) +
# scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)
head(counties)
tmp<-unique(select(filter(Corona_Cases.US_state,Date=="2020-04-17"),c(Lat,Long,Total_confirmed_cases.per100)))
st_as_sf(map("county", plot = FALSE, fill = TRUE))
join::inner_join.sf(Corona_Cases.US_state, counties)
library(sf)
library(sp)
nc <- st_read(system.file("shape/nc.shp", package="sf"))
class(nc)
spdf <- SpatialPointsDataFrame(coords = select(Corona_Cases.US_state,c("Lat","Long")), data = Corona_Cases.US_state,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
head(spdf)
class(spdf)
st_cast(spdf)
filter(Corona_Cases.US_state.summary,Date=="2020-04-20" & Province.State %in% top_states_modified)